The invention relates to a method of determining an intrinsic spectrum f of radiation emitted by an object to be examined, which intrinsic spectrum f is represented by a set of N data points f1 . . . fN and is determined from a measured spectrum h which is represented by a set of M measuring points h1 . . . hM and is measured by means of an analysis apparatus having a given apparatus transfer function G in the form of an Mxc3x97N matrix, which method includes the following steps:
a) forming an approximated intrinsic spectrum g of N data points g1 . . . gN;
b) determining a measure of misfit "khgr"2 between the approximated intrinsic spectrum, convoluted with the apparatus transfer function G, and the measured spectrum;
c) determining the value of a predetermined regularizing function S by inserting the approximated intrinsic spectrum in this function;
d) forming a functional F="khgr"2+xcex1S containing a regularizing constant xcex1;
e) solving the regularizing constant xcex1 from said functional, the N eigenvalues xcex1 . . . xcexN of an Nxc3x97N auxiliary matrix A formed from the apparatus transfer function and the approximated intrinsic spectrum being determined during said solution process;
f) executing a minimizing process on the functional F with the last regularizing constant xcex1 found, using the N data points of the intrinsic spectrum as variables, the N data points thus found constituting a next approximated intrinsic spectrum;
g) repeating the steps b) to f), if necessary, until a predetermined convergence criterion has been satisfied;
h) identifying the approximated intrinsic spectrum then valid as the intrinsic spectrum searched.
The invention also relates to a storage medium which can be read by a computer and is provided with a computer program for carrying out said method, and also to a radiation analysis apparatus which is suitable to carry out the method.
An algorithm for carrying out such a method is known from the article xe2x80x9cBayesian Interpolationxe2x80x9d by David J. C. MacKay, Neural Computation 4, pp. 415-447. Algorithms of this kind are known as xe2x80x9cMaximum Entropy Algorithmsxe2x80x9d.
The cited article, notably chapter 4 thereof, describes how the intrinsic variation of a quantity can be determined from a set of measuring values of the relevant quantity which suffer from noise and other disturbing effects, i.e. the variation of this quantity if all disturbing effects exerted by, for example, the measuring equipment and/or static processes were removed.
A situation of this kind occurs, for example, during the measurement of an intensity spectrum as is done in X-ray diffraction. An object to be examined (a crystalline sample) is then irradiated by X-rays which are emitted again by the sample in a manner which is characteristic of the relevant material. The intensity of such emitted radiation is dependent on the angle at which the radiation is incident on the lattice faces of the crystalline material to be examined. An intensity spectrum of the emitted radiation is measured as a function of the take-off angle by moving an X-ray detector around the sample during the measurement.
In order to achieve a suitable angular resolution of the measurement, a limiting gap with a gap width of the order of magnitude of from 20 to 200 xcexcm is arranged in front of the detector; in the case of a circumscribed circle having a radius of 30 cm this means that a measurement of a spectrum of one revolution yields a number of measuring points N of the order of magnitude of from N=104 to N=105.
As is known from the practice of measurement by means of such apparatus, the spectrum of measuring points is the convolution of the intrinsic spectrum with the apparatus transfer function, possibly increased by noise and contributions by other disturbing effects. The transfer function takes into account the effect of all optical elements in the radiation path from the radiation source to the detector, notably the finite width of the detector slit; furthermore, it is in general also dependent on the location of the measurement (i.e. the magnitude of the take-off angle) so that, as is known for such apparatus, for numerical processing the transfer function takes the form of an Mxc3x97N matrix (where N and M are of the same order of magnitude), i.e. a matrix with from Mxc3x97N=108 to Mxc3x97N=1010 matrix elements.
The determination of the intrinsic spectrum according to the algorithm disclosed in the cited article by MacKay is based on an approximation of the intrinsic spectrum. This approximation may be based on prior theoretical knowledge of the spectrum to be measured, but the measured spectrum consisting of a set of M measuring points may also be taken as the approximated intrinsic spectrum; in order to obtain a set of N points, interpolation can be performed between the M measuring points (if M less than N) or a part of the M measuring points can be ignored (if M greater than N). Subsequently, a measure of misfit "khgr"2 is formed between the approximated intrinsic spectrum, convoluted with the apparatus transfer function, and the measured spectrum. Because of the absence of, for example noise and other statistical functions, this measure of misfit does not have the value zero.
In conformity with the Maximum Entropy Algorithm rule there is formed a functional F="khgr"2+xcex1S in which the regularizing function S is dependent on the N data points of the intrinsic spectrum. The appearance of this regularizing function S is dependent on the nature of the process to be measured; in the case of X-ray diffraction, the appearance of this function may be xcexa3(fi)*log(fi), in which the quantities fi represent the intensities of the measuring points. The regularizing function S includes a factor xcex1 which is referred to as the regularizing constant. According to the Maximum Entropy Algorithm, the functional F must be minimized in dependence on the value of the data points of the intrinsic spectrum. The spectrum of data points at which the functional F is minimum then constitutes the intrinsic spectrum searched. For numerical execution of the minimizing process, however, it is first necessary to determine the value of the regularizing constant xcex1; moreover, an assumption must be made in respect of numerical initial values of the intrinsic spectrum, both in the regularizing function S and in the quantity "khgr"2. As has already been stated, the measured spectrum can often be chosen for the numerical initial values of the intrinsic spectrum.
The cited article describes a process for solving the regularizing constant xcex1 from the functional F. Therein, an Nxc3x97N auxiliary matrix is first determined from the apparatus transfer function and the approximated intrinsic spectrum. The process of forming the auxiliary matrix A is described in chapter 4.3 of the cited article. The set of eigenvalues of this auxiliary matrix is determined. A relation can then be derived between the regularizing constant xcex1 and the set of eigenvalues xcex1 . . . xcexN, so that the regularizing constant xcex1 can be determined from this relation. This process is described notably in chapter 4.4 of the cited article; said relation is found by equating the formulas (4.8) and (4.9) described therein; it then follows that:                               2          ⁢          α          ⁢                      xe2x80x83                    ⁢          S                =                              ∑                          i              =              1                                      i              =              N                                ⁢                                    λ              i                                                      λ                i                            +              α                                                          (        1        )            
(The quantity EWMP stated in the cited article equals the regularizing function S). The value of the regularizing constant xcex1 can then be determined from expression (1) by means of standard solution methods. The value thus found can be inserted in the functional F, after which the minimum value of the functional F is determined in known manner, the N data points of the intrinsic spectrum then being the variables of the functional F. Those values of the data points at which the minim of F occurs constitute a better approximation of the intrinsic spectrum than the values of the initially selected approximated intrinsic spectrum. Depending on the desired precision of the final result, it can then be decided that the approximated intrinsic spectrum thus determined will be the intrinsic spectrum searched; in that case the calculations need not be continued. However, it may also be decided to continue the calculations; the approximated intrinsic spectrum determined thus far then forms the starting point of a next iteration cycle. Such iteration cycles can be repeated until the convergence criterion imposed by the desired precision is satisfied.
The cited article by MacKay does not provide information regarding the dimensions of the auxiliary matrix A. However, it is generally known that the calculation time required to determine the eigenvalues of a matrix increases as the third power of the linear dimension of such a matrix. In the case of matrices such as they may occur, for example in the acquisition of intensity spectra, the value of N or M, so the dimension of the matrix A, may be of the order of magnitude of from 104 to 105. (The quantities M and N are of the same order of magnitude.) The determination of the eigenvalues of a matrix having such dimensions is no longer possible in practical circumstances. Notably for analytic equipment as used for standard laboratory analysis, it is desirable to execute the associated calculations with a customary personal computer with a calculation time which is of the same order of magnitude as the time required for executing the associated measurement.
It is an object of the invention to provide a method of the kind set forth in which the determination of the eigenvalues of matrices of large dimensions is significantly faster than in the standard prior art methods.
To this end, the method according to the invention is characterized in that for the determination of the N eigenvalues xcex1 . . . xcexN of the Nxc3x97N auxiliary matrix A, the auxiliary matrix is subdivided into a number of L partial auxiliary matrices Pj (j=1 . . . L) which are situated around the diagonal of the auxiliary matrix A,
and that the eigenvalues of each partial auxiliary matrix Pj are determined separately,
the set of eigenvalues xcex1 . . . xcexN of the Nxc3x97N auxiliary matrix A consisting of the set of all eigenvalues of the partial auxiliary matrices Pj.
The invention is based on the recognition of the fact that the apparatus transfer function G has the form of a Mxc3x97N matrix (where M is the number of measuring points) for many spectroscopic measurements, but that the numbers in this matrix are much greater in the direct vicinity of the diagonal than the numbers which are situated further away. This is due to the fact that the width of the detector slit of such a spectroscopic apparatus is generally very small relative to the total measuring trajectory, so that a point of the intrinsic spectrum which is situated far from the slit position has only a (very) slight effect on the intensity measured in the relevant slit position. Thus, only a narrow band of numbers around the diagonal is of importance in the rows of the matrix.
The formation of the Nxc3x97N auxiliary matrix A from the apparatus transfer function and the approximated intrinsic spectrum is performed in such a manner that the appearance of this auxiliary matrix has a structure which is comparable to that of the apparatus transfer function G, so that A also has a band structure.
The invention is also based on the recognition of the fact that the determination of the eigenvalues of such an auxiliary matrix A can be performed by replacing said auxiliary matrix by a comparatively large number of much smaller partial auxiliary matrices P which succeed one another along the diagonal of the auxiliary matrix and whose diagonals lie on that of the auxiliary matrix. The dimensions of such partial auxiliary matrices are chosen to be such that the elements in each row of said partial auxiliary matrices contain the large values of the corresponding row of the auxiliary matrix and that, consequently, the elements of the relevant row which lie outside the partial auxiliary matrix have a value such that it can be ignored for all practical purposes. The determination of the set of eigenvalues of such a (much) smaller partial auxiliary matrix is much easier than that of the large auxiliary matrix. The gain in respect of calculation time can be illustrated on the basis of the following numerical example: when the auxiliary matrix A is subdivided into 100 partial auxiliary matrices P, the calculation time for determining the eigenvalues of one partial auxiliary matrix becomes 106 times shorter, so 104 times shorter for all 100 partial auxiliary matrices. A substantial gain is thus realized in respect of the calculation time required to determine the eigenvalues.
By performing an experimental check on the basis of a known spectrum (or by once applying a conventional, time-consuming diagonalizing algorithm to the auxiliary matrix), it can be checked whether the precision thus achieved in determining the eigenvalues of the auxiliary matrix A suffices for practical purposes. Such adequate precision can also be demonstrated theoretically.
In a version of the invention the eigenvalues of the partial auxiliary matrices are determined by determination of the Fourier transform of an arbitrary row of each of said partial auxiliary matrices.
The insight on which this aspect of the invention is based consists in that said narrow band of elements along the diagonal of the auxiliary matrix varies only slowly in value, i.e. even though the numbers in said band of numbers have been shifted one location between two successive rows, these numbers exhibit only small differences per row. This small difference per row is due to the fact that, because of the small width of the detector slit, the transfer from the emissive sample to the detector varies only gradually as a function of the angular position of the detector. A matrix exhibiting both these phenomena (i.e. a band structure and a small difference between the elements of neighboring rows along the diagonal) resembles a so-called Toeplitz matrix. As is known, a Toeplitz matrix is a matrix whose first row consists of arbitrary numbers. The subsequent row is obtained by shifting all elements of the previous row through one location in the row; the void then arising at one end is filled with an arbitrary element and the element at the other end is deleted. (When the void arising at one end is filled with the element deleted at the other end, a so-called xe2x80x9ccirculant continuationxe2x80x9d is concerned). This phenomena is repeated for all subsequent rows of the matrix.
The auxiliary matrix is distinct from a Toeplitz matrix in that a given variation of the value of the elements occurs along the matrix diagonal in the auxiliary matrix. As has already been described, the determination of the eigenvalues of such an auxiliary matrix is performed by replacing said auxiliary matrix by a comparatively large number of much smaller partial auxiliary matrices which succeed one another along the diagonal of the auxiliary matrix and whose diagonals lie on that of the auxiliary matrix. The dimensions of said partial auxiliary matrices must be chosen to be such that the variation of the values of the elements along their diagonal is negligibly small for all practical purposes. Because such partial auxiliary matrices approximate the Toeplitz structure much better than the auxiliary matrix itself, a theorem which holds for Toeplitz matrices can be applied to such partial auxiliary matrices; according to theorem the eigenvalues of a Toeplitz matrix with a circulant continuation exactly equal the coefficients of the Fourier transform of an arbitrary row of the relevant matrix. Even though the partial auxiliary matrices need not be subjected to a circulant continuation in practice, the eigenvalues obtained by Fourier transformation constitute a suitable approximation of the exact eigenvalues of the partial auxiliary matrices. To those skilled in the art it will be generally known that the calculation time for determining the eigenvalues by means of Fourier transformation corresponds to N*log(N) (in which N is the dimension of the matrix) instead of N3 as in the case of conventional determination of eigenvalues. The determination of the eigenvalues of the partial auxiliary matrices is thus reduced to the determination of the Fourier coefficients which is far less intensive in respect of calculation work, i.e. to the determination of the discrete Fourier transform of each time one matrix row of each of the partial matrices.
Ignoring the values in the auxiliary matrix which are not included in a partial auxiliary matrix, introduces a first error in the determination of the eigenvalues. Moreover, the determination of the eigenvalues of the partial auxiliary matrices by Fourier transformation (actually by equating the partial auxiliary matrices with the circulant continuation of a Toeplitz matrix) introduces a second error. A further advantage achieved by this version of the invention consists in that the first and the second error oppose one another, so that the resultant error is smaller than each of said errors individually.
According to a further version of the invention, the eigenvalues of the partial auxiliary matrices are determined by determining the Fourier transform of the mean of at least two arbitrary rows of each of said partial auxiliary matrices.
Because of this step, the actual eigenvalues are approximated even better since a given variation of the values of elements of the different rows is compensated because a mean row is thus used.
According to a further version of the invention, an amount which is arbitrarily distributed among the measuring points is added to the measuring points of the set of M measuring points during the execution of the minimizing process.
During the determination of the minimum of the functional F it may occur that, after a given iteration step in the minimizing process, the variation of this functional around the minimum is such that during a next iteration step either an extremely small further approximation of the minimum is realized (in the case of a very flat variation of the functional in said range, or that during a next iteration step the minimum is passed and the process is returned to the original iteration point again by the subsequent iteration step (in the case of a very acute variation of the functional in said range). In both cases the approximation of the minimum is very slow or even non-existent. Such undesirable phenomena can be avoided by the addition of artificial noise which decreases to zero as the minimizing process progresses.